# Loading the csv files from dropbox
hamilton_final_price <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_price_old.csv", quote="\"")
hamilton_final_price_for7 <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_price.csv", quote="\"")
hamilton_final_CornS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_CornS.csv", quote="\"")
hamilton_final_CornF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_CornF.csv", quote="\"")
hamilton_final_OatsS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_OatsS.csv", quote="\"")
hamilton_final_OatsF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_OatsF.csv", quote="\"")
hamilton_final_RyeS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_RyeS.csv", quote="\"")
hamilton_final_RyeF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_RyeF.csv", quote="\"")
hamilton_final_LardS <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_LardS.csv", quote="\"")
hamilton_final_LardF <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_LardF.csv", quote="\"")
hamilton_final_TBill <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_TBill.csv", quote="\"")
hamilton_final_RyeS <- hamilton_final_RyeS[-58]
hamilton_final_RyeF <- hamilton_final_RyeF[-58]
hamilton_final_SepDummy <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_SepDummy.csv", quote="\"")
hamilton_final_MayDummy <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_MayDummy.csv", quote="\"")
# converting them into time series objects
Sep.ts = ts(data=hamilton_final_SepDummy, frequency = 3, start=c(1921,1), end=c(1939,1)) 
May.ts = ts(data=hamilton_final_MayDummy, frequency = 3, start=c(1921,1), end=c(1939,1)) 
# converting them into time series objects
Price.ts = ts(data=hamilton_final_price, frequency = 3, start=c(1920,2), end=c(1939,2)) 
Price_for7.ts= ts(data=hamilton_final_price_for7, frequency = 3, start=c(1920,1), end=c(1939,2)) 

CornS.ts = ts(data=hamilton_final_CornS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
CornF.ts = ts(data=hamilton_final_CornF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
OatsS.ts = ts(data=hamilton_final_OatsS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
OatsF.ts = ts(data=hamilton_final_OatsF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
RyeS.ts = ts(data=hamilton_final_RyeS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
RyeF.ts = ts(data=hamilton_final_RyeF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
LardS.ts = ts(data=hamilton_final_LardS, frequency = 3, start=c(1920,2), end=c(1939,2)) 
LardF.ts = ts(data=hamilton_final_LardF, frequency = 3, start=c(1920,2), end=c(1939,2)) 
TBill.ts = ts(data=hamilton_final_TBill, frequency = 3, start=c(1920,2), end=c(1939,2)) 

# To extract F p-values 
lmp <- function (modelobject) {
  if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
  f <- summary(modelobject)$fstatistic
  p <- pf(f[1],f[2],f[3],lower.tail=F)
  attributes(p) <- NULL
  return(p)
}
# To extract F values
fvlu <- function (modelobject) {
  f <- summary(modelobject)$fstatistic[["value"]]
  return(f)
}
fvlu(Corn1.lm)



# Hamilton cares about the (forward) lagged differences between spot and futures, so below I examine these
FwdCornS.ts <- lag(CornS.ts, k=1)
FwdOatsS.ts <- lag(OatsS.ts, k=1)
FwdRyeS.ts <- lag(RyeS.ts, k=1)
FwdLardS.ts <- lag(LardS.ts, k=1)
# From now on, we are in t+1 space

CornSpotPred.ts <- 100 *(log(FwdCornS.ts)-log(CornF.ts))
OatsSpotPred.ts <- 100 *(log(FwdOatsS.ts)-log(OatsF.ts))
RyeSpotPred.ts <- 100 *(log(FwdRyeS.ts)-log(RyeF.ts))
LardSpotPred.ts <- 100 *(log(FwdLardS.ts)-log(LardF.ts))

# Corn Errors v and v(t-1) (which is really like one and two period lags on V (t+1))

LagErrorCorn1.ts <- lag(CornSpotPred.ts,k=-1)
LagErrorCorn2.ts <- lag(LagErrorCorn1.ts,k=-1)
LagErrorCorn1Adj.ts <- LagErrorCorn1.ts[-c(1,56)]
LagErrorCorn2Adj.ts <- LagErrorCorn2.ts[-c(55,56)]
LagErrorCorn1Final.ts = ts(LagErrorCorn1Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
LagErrorCorn2Final.ts = ts(LagErrorCorn2Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
CornSpotPred1.ts <- CornSpotPred.ts[-c(1,2)]
# Check that the data are aligned
CornLags.ts <- cbind(CornSpotPred1.ts, LagErrorCorn1Final.ts, LagErrorCorn2Final.ts)
head(CornLags.ts)
Corn1.lm <- lm(CornSpotPred1.ts ~ LagErrorCorn1Final.ts + LagErrorCorn2Final.ts)
summary(Corn1.lm)

TBillMatchedUnadj.ts <- TBill.ts[-58]
CornSpotPred.ts = ts(CornSpotPred.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
TBillMatched.ts = ts(TBillMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Corn2.lm <- lm(CornSpotPred.ts ~ TBillMatched.ts)

# summary(Corn2.lm)
# CornTBill.ts <- cbind(CornSpotPred.ts, TBillMatched.ts)
# head(CornTBill.ts)

CornFMatchedUnadj.ts <- (log(CornF.ts)*100)[-58]
CornFMatched.ts = ts(CornFMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Corn3.lm <- lm(CornSpotPred.ts ~ CornFMatched.ts)
summary(Corn3.lm)
CornCornS.ts <- cbind(CornSpotPred.ts, CornFMatched.ts)
head(CornCornS.ts)

CornSMatchedUnadj.ts <- (log(CornS.ts)*100)[-58]
CornSMatched.ts = ts(CornSMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Corn4.lm <- lm(CornSpotPred.ts ~ CornSMatched.ts)
summary(Corn4.lm)
CornCornS.ts <- cbind(CornSpotPred.ts, CornSMatched.ts)
head(CornCornS.ts)

SeasC1 <- cycle(CornSpotPred.ts)
Corn5.lm <- lm(CornSpotPred.ts ~ factor(SeasC1))
summary(Corn5.lm)

# Now it gets complicated with all of the above variables 
SeasC2 <- cycle(CornSpotPred1.ts)
TBillMatched1.ts <- TBillMatched.ts[-c(1,2)]
CornFMatched1.ts <- CornFMatched.ts[-c(1,2)]
CornSMatched1.ts <- CornSMatched.ts[-c(1,2)]
Corn6.lm <- lm(CornSpotPred1.ts ~ LagErrorCorn1Final.ts + LagErrorCorn2Final.ts + TBillMatched1.ts + CornFMatched1.ts + CornSMatched1.ts + Sep.ts + May.ts)
summary(Corn6.lm)

# Leave the final 7th set of regressions until the end 
# Now Oats

LagErrorOats1.ts <- lag(OatsSpotPred.ts,k=-1)
LagErrorOats2.ts <- lag(LagErrorOats1.ts,k=-1)
LagErrorOats1Adj.ts <- LagErrorOats1.ts[-c(1,56)]
LagErrorOats2Adj.ts <- LagErrorOats2.ts[-c(55,56)]
LagErrorOats1Final.ts = ts(LagErrorOats1Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
LagErrorOats2Final.ts = ts(LagErrorOats2Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
OatsSpotPred1.ts <- OatsSpotPred.ts[-c(1,2)]
# Check that the data are aligned
# OatsLags.ts <- cbind(OatsSpotPred1.ts, LagErrorOats1Final.ts, LagErrorOats2Final.ts)

Oats1.lm <- lm(OatsSpotPred1.ts ~ LagErrorOats1Final.ts + LagErrorOats2Final.ts)
summary(Oats1.lm)

TBillMatchedUnadj.ts <- TBill.ts[-58]
OatsSpotPred.ts = ts(OatsSpotPred.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
TBillMatched.ts = ts(TBillMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Oats2.lm <- lm(OatsSpotPred.ts ~ TBillMatched.ts)
summary(Oats2.lm)
# OatsTBill.ts <- cbind(OatsSpotPred.ts, TBillMatched.ts)


OatsFMatchedUnadj.ts <- (log(OatsF.ts)*100)[-58]
OatsFMatched.ts = ts(OatsFMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Oats3.lm <- lm(OatsSpotPred.ts ~ OatsFMatched.ts)
summary(Oats3.lm)
# OatsOatsF.ts <- cbind(OatsSpotPred.ts, OatsFMatched.ts)


OatsSMatchedUnadj.ts <- (log(OatsS.ts)*100)[-58]
OatsSMatched.ts = ts(OatsSMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Oats4.lm <- lm(OatsSpotPred.ts ~ OatsSMatched.ts)
summary(Oats4.lm)
# OatsOatsS.ts <- cbind(OatsSpotPred.ts, OatsSMatched.ts)


SeasO1 <- cycle(OatsSpotPred.ts)
Oats5.lm <- lm(OatsSpotPred.ts ~ factor(SeasO1))
summary(Oats5.lm)

# Now all of the above variables 
TBillMatched1.ts <- TBillMatched.ts[-c(1,2)]
OatsFMatched1.ts <- OatsFMatched.ts[-c(1,2)]
OatsSMatched1.ts <- OatsSMatched.ts[-c(1,2)]
Oats6.lm <- lm(OatsSpotPred1.ts ~ LagErrorOats1Final.ts + LagErrorOats2Final.ts + TBillMatched1.ts + OatsFMatched1.ts + OatsSMatched1.ts + Sep.ts + May.ts)
summary(Oats6.lm)

# Now Rye

LagErrorRye1.ts <- lag(RyeSpotPred.ts,k=-1)
LagErrorRye2.ts <- lag(LagErrorRye1.ts,k=-1)
LagErrorRye1Adj.ts <- LagErrorRye1.ts[-c(1,56)]
LagErrorRye2Adj.ts <- LagErrorRye2.ts[-c(55,56)]
LagErrorRye1Final.ts = ts(LagErrorRye1Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
LagErrorRye2Final.ts = ts(LagErrorRye2Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
RyeSpotPred1.ts <- RyeSpotPred.ts[-c(1,2)]
Rye1.lm <- lm(RyeSpotPred1.ts ~ LagErrorRye1Final.ts + LagErrorRye2Final.ts)

TBillMatchedUnadj.ts <- TBill.ts[-58]
RyeSpotPred.ts = ts(RyeSpotPred.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
TBillMatched.ts = ts(TBillMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Rye2.lm <- lm(RyeSpotPred.ts ~ TBillMatched.ts)


RyeFMatchedUnadj.ts <- (log(RyeF.ts)*100)[-58]
RyeFMatched.ts = ts(RyeFMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Rye3.lm <- lm(RyeSpotPred.ts ~ RyeFMatched.ts)


RyeSMatchedUnadj.ts <- (log(RyeS.ts)*100)[-58]
RyeSMatched.ts = ts(RyeSMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Rye4.lm <- lm(RyeSpotPred.ts ~ RyeSMatched.ts)


SeasO1 <- cycle(RyeSpotPred.ts)
Rye5.lm <- lm(RyeSpotPred.ts ~ factor(SeasO1))


TBillMatched1.ts <- TBillMatched.ts[-c(1,2)]
RyeFMatched1.ts <- RyeFMatched.ts[-c(1,2)]
RyeSMatched1.ts <- RyeSMatched.ts[-c(1,2)]
Rye6.lm <- lm(RyeSpotPred1.ts ~ LagErrorRye1Final.ts + LagErrorRye2Final.ts + TBillMatched1.ts + RyeFMatched1.ts + RyeSMatched1.ts + Sep.ts + May.ts)


# Now Lard

LagErrorLard1.ts <- lag(LardSpotPred.ts,k=-1)
LagErrorLard2.ts <- lag(LagErrorLard1.ts,k=-1)
LagErrorLard1Adj.ts <- LagErrorLard1.ts[-c(1,56)]
LagErrorLard2Adj.ts <- LagErrorLard2.ts[-c(55,56)]
LagErrorLard1Final.ts = ts(LagErrorLard1Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
LagErrorLard2Final.ts = ts(LagErrorLard2Adj.ts, frequency = 3, start=c(1921,1), end=c(1939,1)) 
LardSpotPred1.ts <- LardSpotPred.ts[-c(1,2)]
Lard1.lm <- lm(LardSpotPred1.ts ~ LagErrorLard1Final.ts + LagErrorLard2Final.ts)


TBillMatchedUnadj.ts <- TBill.ts[-58]
LardSpotPred.ts = ts(LardSpotPred.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
TBillMatched.ts = ts(TBillMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Lard2.lm <- lm(LardSpotPred.ts ~ TBillMatched.ts)


LardFMatchedUnadj.ts <- (log(LardF.ts)*100)[-58]
LardFMatched.ts = ts(LardFMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Lard3.lm <- lm(LardSpotPred.ts ~ LardFMatched.ts)


LardSMatchedUnadj.ts <- as.ts((log(LardS.ts)*100)[-58])
LardSMatched.ts = ts(LardSMatchedUnadj.ts, frequency = 3, start=c(1920,2), end=c(1939,1)) 
Lard4.lm <- lm(LardSpotPred.ts ~ LardSMatched.ts)

SeasL1 <- cycle(LardSpotPred.ts)
Lard5.lm <- lm(LardSpotPred.ts ~ factor(SeasL1))

TBillMatched1.ts <- TBillMatched.ts[-c(1,2)]
LardFMatched1.ts <- LardFMatched.ts[-c(1,2)]
LardSMatched1.ts <- LardSMatched.ts[-c(1,2)]
Lard6.lm <- lm(LardSpotPred1.ts ~ LagErrorLard1Final.ts + LagErrorLard2Final.ts + TBillMatched1.ts + LardFMatched1.ts + LardSMatched1.ts + Sep.ts + May.ts)

# Now the 7th set of regressions in table 2
# Special matching - use may1 and Sept dummy files
hamilton_final_Sep1Dummy <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_SepDummy1.csv", quote="\"")
hamilton_final_May1Dummy <- read.table("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Data/OriginalData/hamilton_final_MayDummy1.csv", quote="\"")
# converting them into time series objects
Sep1.ts = ts(data=hamilton_final_Sep1Dummy, frequency = 3, start=c(1920,2), end=c(1939,1)) 
May1.ts = ts(data=hamilton_final_May1Dummy, frequency = 3, start=c(1920,2), end=c(1939,1)) 
# I have fixed the price data to add one more data point, so I need to match the data again
PriceMatched.ts <- log(Price_for7.ts[-c(0,59)])*100
LagPrice.ts <-lag(Price_for7.ts,-1)
LagPriceMatched.ts <- log(LagPrice.ts[-c(0,58,59)])*100
PriceMatched.ts <- PriceMatched.ts[-c(1)]
# The next eight lines are not technically necessary
CornSpotPred2.ts <- CornSpotPred.ts
OatsSpotPred2.ts <- OatsSpotPred.ts
RyeSpotPred2.ts <- RyeSpotPred.ts
LardSpotPred2.ts <- LardSpotPred.ts
CornSMatched2.ts <- CornSMatched.ts
OatsSMatched2.ts <- OatsSMatched.ts
RyeSMatched2.ts <- RyeSMatched.ts
LardSMatched2.ts <- LardSMatched.ts

Corn7.lm <- lm(CornSpotPred2.ts ~ PriceMatched.ts + LagPriceMatched.ts + CornSMatched2.ts + OatsSMatched2.ts +  RyeSMatched2.ts + LardSMatched2.ts + Sep1.ts + May1.ts)
summary(Corn7.lm)

Oats7.lm <- lm(OatsSpotPred2.ts ~ PriceMatched.ts + LagPriceMatched.ts + CornSMatched2.ts + OatsSMatched2.ts +  RyeSMatched2.ts + LardSMatched2.ts + Sep1.ts + May1.ts)
summary(Oats7.lm)

Rye7.lm <- lm(RyeSpotPred2.ts ~ PriceMatched.ts + LagPriceMatched.ts + CornSMatched2.ts + OatsSMatched2.ts +  RyeSMatched2.ts + LardSMatched2.ts + Sep1.ts + May1.ts)
summary(Rye7.lm)

Lard7.lm <- lm(LardSpotPred2.ts ~ PriceMatched.ts + LagPriceMatched.ts + CornSMatched2.ts + OatsSMatched2.ts +  RyeSMatched2.ts + LardSMatched2.ts + Sep1.ts + May1.ts)
summary(Lard7.lm)

sink("~/Dropbox/Replication Workshop/Student Papers/Saleuddin/Analysis/R/table2.out")
# Replicate Table 2
"Table 2 - Market Efficiency Tests"

Table1 <- matrix (c(summary(Corn1.lm)$fstatistic[["dendf"]], summary(Corn1.lm)$fstatistic[["numdf"]],fvlu(Corn1.lm),lmp(Corn1.lm),fvlu(Oats1.lm),lmp(Oats1.lm), fvlu(Rye1.lm), lmp(Rye1.lm), fvlu(Lard1.lm),lmp(Lard1.lm)
                    , summary(Corn2.lm)$fstatistic[["dendf"]], summary(Corn2.lm)$fstatistic[["numdf"]],fvlu(Corn2.lm),lmp(Corn2.lm),fvlu(Oats2.lm),lmp(Oats2.lm), fvlu(Rye2.lm), lmp(Rye2.lm), fvlu(Lard2.lm),lmp(Lard2.lm)
                    , summary(Corn3.lm)$fstatistic[["dendf"]], summary(Corn3.lm)$fstatistic[["numdf"]],fvlu(Corn3.lm),lmp(Corn3.lm),fvlu(Oats3.lm),lmp(Oats3.lm), fvlu(Rye3.lm), lmp(Rye3.lm), fvlu(Lard3.lm),lmp(Lard3.lm)
                    , summary(Corn4.lm)$fstatistic[["dendf"]], summary(Corn4.lm)$fstatistic[["numdf"]],fvlu(Corn4.lm),lmp(Corn4.lm),fvlu(Oats4.lm),lmp(Oats4.lm), fvlu(Rye4.lm), lmp(Rye4.lm), fvlu(Lard4.lm),lmp(Lard4.lm)
                    , summary(Corn5.lm)$fstatistic[["dendf"]], summary(Corn5.lm)$fstatistic[["numdf"]],fvlu(Corn5.lm),lmp(Corn5.lm),fvlu(Oats5.lm),lmp(Oats5.lm), fvlu(Rye5.lm), lmp(Rye5.lm), fvlu(Lard5.lm),lmp(Lard5.lm)
                    , summary(Corn6.lm)$fstatistic[["dendf"]], summary(Corn6.lm)$fstatistic[["numdf"]],fvlu(Corn6.lm),lmp(Corn6.lm),fvlu(Oats6.lm),lmp(Oats6.lm), fvlu(Rye6.lm), lmp(Rye6.lm), fvlu(Lard6.lm),lmp(Lard6.lm)
                    , summary(Corn7.lm)$fstatistic[["dendf"]], summary(Corn7.lm)$fstatistic[["numdf"]],fvlu(Corn7.lm),lmp(Corn7.lm),fvlu(Oats7.lm),lmp(Oats7.lm), fvlu(Rye7.lm), lmp(Rye7.lm), fvlu(Lard7.lm),lmp(Lard7.lm)
                    ), ncol=10, byrow=TRUE)
colnames(Table1) <- c("F dendf", "F numdf", "Corn F", "Corn P", "Oats F", "Oats P", "Rye F", "Rye P", "Lard F", "Lard P")
rownames(Table1) <- c("Lagged errors", "Interest rate", "Futures price", "Spot price", "Seasonal dummies", "All above","All VAR")
# rownames(Table1) <- c("Lagged errors t and t-1","Interest rate (t)", "Futures price (t)", "Spot price(t)","Seasonal dummies d1 and d2","All the above", "All VAR variables")
Table1 <- as.table(Table1)
Table1

sink()
